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ARC-LENGTH  CONTINUATION  AND  MULTI-GRID  TECHNIQUES 
FOR  NONLINEAR  ELLIPTIC  EIGENVALUE  PROBLEMS" 

TONY  F.  C.  CHAN'*  and  H.  B.  KELLER} 

Abstract.  We  investigate  multi-grid  methods  for  solving  linear  systems  arising  from  arc-length  continu¬ 
ation  technique*  applied  to  nonlinear  elliptic  eigenvalue  problems.  We  find  thai  the  usual  multi-grid  method* 
diverge  in  the  neighborhood  of  singular  points  of  the  solution  branches.  As  a  result,  the  continuation 
method  is  unable  to  continue  past  a  limit  point  in  the  Bratu  problem.  This  divergence  is  analyzed  and  a 
modified  multi-grid  algorithm  has  been  devised  based  on  this  analysis  In  principle,  this  new  multi-grid 
algorithm  converges  for  elliptic  systems  arbitrarily  close  to  singularity  and  has  been  used  successfully  in 
conjunction  with  arc-length  continuation  procedures  on  the  model  problem.  In  the  worst  situation,  both 
the  storage  and  the  computational  work  are  only  about  a  factor  of  two  more  than  the  unmodified  multi-grid 
methods. 
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1.  Introduction.  Many  problems  of  computational  interest  can  be  formulated  as 
(1.1)  G(u.A)-0, 

where  u  represents  the  “solution"  <i.e.,  flow  field,  displacements,  etc.)  and  A  is  a  real 
physical  parameter  (i.e.,  Reynold's  number,  load,  etc.)  It  is  required  to  find  the  solution 
for  some  A  -intervals,  that  is.  a  path  of  solutions,  [i#(A ),  A  ].  In  this  paper,  we  use  a 
class  of  continuation  based  on  parametrizing  the  solution  branches  by  arc-length,  say 
[«(5),  A  (5)].  A  main  advantage  of  these  arc-length  continuation  methods  is  that  most 
singular  points  on  the  solution  branches  can  be  handled  without  much  difficulty. 
Equations  of  the  form  1 1.1 1  are  called  nonlinear  elliptic  eigenvalue  problems  if  the 
operator  G  with  A  fixed  is  an  elliptic  differential  operator  [2].  For  nonlinear  elliptic 
eigenvalue  problems,  a  major  portion  of  the  computational  work  in  the  arc-length 
continuation  methods  is  spent  in  solving  large  linear  elliptic  systems.  In  this  paper, 
we  investigate  the  use  of  multi-grid  [4]  methods  for  solving  these  linear  systems.  It 
turns  out  that  a  straightforward  implementation  of  the  multi-grid  methods  fails  in  the 


neighborhood  of  the  singular  points  and  this  usually  prevents  continuation  past  limit 
points.  This  failure  is  analyzed  and  a  modified  multi-grid  method  based  on  this  analysis 
is  devised.  Even  for  very  singular  systems,  the  new  multi-grid  algorithm  performs 
satisfactorily  and  never  requires  more  than  about  twice  the  storage  and  computational 
work  as  the  unmodified  algorithm. 

The  arc-length  continuation  methods  will  be  described  in  $  2  and  the  multi-grid 
methods  in  S  3.  In  $  4,  computational  results  for  a  model  problem  are  presented, 
together  with  a  description  of  the  difficulties  encountered  by  the  multi-grid  method 
near  a  limit  point.  The  behavior  of  the  multi-grid  method  near  singular  points  will  be 
analyzed  in  $  S.  The  modified  multi-grid  algorithms  designed  to  overcome  these 
difficulties  are  described  in  6  6.  The  paper  ends  with  a  summary  in  $  7. 
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2.  Newton’s  method  and  continuation  techniques.  In  this  paper  we  are  concerned 
with  methods  for  computing  a  family  or  path  of  solutions  of  (1.1).  The  methods  we 
employ  will  be  based  on  some  version  of  Newton's  method. 

2.1.  Newton’s  method.  Given  a  value  of  A  and  an  initial  guess  u°  for  the  solution 
i/(A),  we  perform  the  following  steps  repeatedly  until  ||Su'||<  e  is  satisfied: 

(2.1)  GL«u'--G(«'.A). 

(2.2)  u'*'  «  u1  +  Su'. 

In  the  above,  subscripts  denote  partial  derivatives  and  so  G„  denotes  the  Jacobian  of 
the  operator  G  (with  respect  to  u ).  This  procedure  will  generally  converge  quadratically 
when  it  does  converge.  However,  as  is  well  known,  in  many  instances  it  will  fail  to 
converge  when  the  initial  guess  is  not  “dose"  to  the  true  solution. 

2.2.  Natural  continuation.  A  plausible  procedure  for  overcoming  this  conver¬ 
gence  difficulty  and  also  for  determining  the  dependence  of  u  on  A  is  to  start  at  a 
known  solution  ( w0 ,  A0)  on  the  solution  curve  and  use  it  as  initial  guess  for  a  Newton- 
type  iteration  to  find  the  solution  for  a  neighboring  point  on  the  solution  curve  with 
A  close  to  Ac-  The  procedure  is  then  repeated.  We  can  improve  on  this  by  computing 
the  derivative,  u*,  at  a  known  solution  and  use  it  to  get  a  better  initial  guess  for  the 
next  value  of  A  in  a  predictor -corrector  fashion.  We  call  this  a  natural  continuation 
procedure  because  it  corresponds  to  parametrizing  the  solution  curve  by  A,  the  naturally 
occurring  parameter.  A  specific  form  of  this  is  the  more  or  less  well-known 

Euler-Newton  continuation  procedure.  Given  a  known  solution  (u0,  A0),  we  com¬ 
pute  the  solutions  at  nearby  values  of  X  as  follows : 

1.  First  compute  the  derivative  u»  at  <u0,  A0)  from 

(2.3)  GHuk  -  -Ga. 

2.  Perform  an  Euler  predictor  step  : 

(2.4)  u°  =  u0  +  «a  (A -A0). 

3.  Use  u°  as  initial  guess  in  Newton 's  method, 

(2.5)  GUm'*'  ~  «f)  “  -G(u‘,  A ), 
until  convergence. 

4.  Use  («( A ),  A )  as  the  new  (uo,  X0)  and  go  back  to  Step  1. 

Note  that  the  computation  of  the  derivative  uk  does  not  cause  much  computational 
overhead  because  we  usually  have  the  factorization  of  the  Jacobian  G„  computed 
already  in  the  Newton  step.  Using  such  a  predictor-corrector  method  will  often  allow 
us  to  take  a  much  bigger  step  in  A  and  thus  reduce  the  overall  cost  of  determining 
the  dependence  of  u  on  A. 

Unfortunately,  this  procedure  needs  some  modification  in  order  to  handle  general 
nonlinear  systems  because  of  the  possibility  of  existence  of  nonunique  solutions.  The 
nonuniqueness  usually  manifests  itself  in  the  form  of  existence  of  "singular”  points 
where  the  Jacobian  G„  is  singular  (see  Fig.  2.1).  Points  such  as  point  A  in  Fig.  2.1 
are  called  limit  points  (or  turning  points)  and  points  such  as  point  B  are  called 
bifurcation  points.  These  singular  points  are  further  characterized  by  the  conditions 
that  G*<£  Range  (G„>  at  a  limit  point  and  that  GA  €  Range  (G„)  at  a  bifurcation  point 
[12]. 

The  difficulties  that  a  natural  continuation  procedure  will  encounter  at  singular 
points  are  threefold.  First  of  all,  since  G.  is  singular  at  these  points,  Newton's  method 


Fig.  2.1.  A  typical  bifurcation  diagram. 
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will  at  best  be  linearly  convergent,  making  it  much  more  costly  to  compute  the  solution. 
Moreover,  near  a  limit  point,  there  may  not  exist  a  solution  for  a  given  value  of  A 
(see  Fig.  2.2)  and  hence  the  iterations  must  fail  to  converge.  Lastly,  we  need  some 
mechanism  for  switching  branches  at  a  bifurcation  point. 

2J.  Arc-length  continuation.  In  the  pseudo  arc-length  continuation  approach 
[12],  these  difficulties  are  overcome  by  not  parametrizing  the  solution  u  by  A.  Instead, 
we  parametrize  the  solution  branches  using  an  arc-length  parameter  s,  and  specify 
how  far  along  the  current  solution  branch  we  want  to  march. 

To  be  more  specific,  we  let  s  be  the  arc-length  parameter,  and  treat  u(s)  and 
Air)  as  functions  of  s.  We  can  compute  the  "tangent"  Ms0)]  (where  the  dots 
denote  differentiation  w  ith  respect  to  s)  of  a  known  solution  at  s  *  r0  from  the  following 
two  equations. 

(2.6i  G„iio  *•  A0G*  =  0, 

12.7)  iW::*|A0iJ-l=0. 

Equation  < 2.6 1  is  obtained  from  differentiating  G(u,  A )  *  0  with  respect  to  s  and  (2.7) 
imposes  the  arc-length  condition.  We  could  theoretically  generate  the  solution  curve 
by  integrating  the  initial  value  problem  obtained  by  solving  (2.6),  (2.7)  for  tits )  and 
A(s  i.  Although  this  process  is  subject  to  the  usual  instabilities  inherent  in  solving  initial 
value  problems  approximately,  it  can  be  an  extremely  effective  procedure.  Indeed  our 
pseudo  arc-length  continuation  procedure  can  be  viewed  as  a  method  for  stabilizing 
Euler  integration  of  l2.6).  (2.7). 


solution  curve  on  which 
C<U(S!.  AtS  II  “0 


Fig. 


Pseudo  arc-leng'n  continuation. 
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In  the  pseudo  aic-length  continuation  procedure,  we  advance  from  s0tos  along 
the  tangent  to  the  solution  brand)  and  require  the  new  solution  u  (r)  and  A  (s)  to  satisfy 

(2.8)  N(i/(*),A(5))«iio(«(*)“*Uo))+Ao(A(j)-A(*o))-(»-*o)“0. 

In  addition  we  require,  of  course, 

(2.9)  G(«(s).  A(s))-0. 

Equation  (2.8)  is  the  linearization  of  (2.7),  and  as  indicated  forces  the  new  solution 
to  lie  on  a  hyperplane  perpendicular  to  the  tangent  vector  to  the  solution  curve  at  s0 
and  at  a  distance  (s-io)  from  it.  Equation  (2.9)  requires  ii(s)  and  A(s)  to  lie  on  the 
true  solution  curve  (Fig.  2.3).  We  now  solve  the  coupled  system  (2.8)  and  (2.9)  for 
u(s)  and  A(r),  given  the  step  size  (s-io)  (efficient  strategies  for  choosing  the  step 
size  are  discussed  in  [23]).  We  use  Newton's  method,  in  which  case  we  have  to  solve 
the  following  linear  system  at  each  iteration: 


It  can  be  shown  that  at  limit  points,  where  G„  is  singular  and  G*/ Range  (GJ, 
the  linear  system  in  (2.10)  is  nonsingular  (see  [12])  and  therefore  Newton's  method 
for  the  coupled  system  (2.8)  and  (2.9)  is  well  defined.  Hence  limit  points  present  no 
problem  and  even  quadratic  convergence  is  achievable. 

At  bifurcation  points,  where  G„  is  singular  and  G*  c  Range  (G„),  things  are  more 
complicated.  In  the  simplest  case  of  only  one  brand)  bifurcating  from  the  main  branch 
(simple  bifurcation),  an  additional  higher  order  condition  involving  G*.  G«*  and  G*» 
has  to  be  satisfied.  It  can  be  shown  [12]  that  this  condition,  together  with  (2.6)  and 
(2.7)  and  the  left  and  right  null  vectors  of  Cm  enable  two  solutions  for  («i0,  Ao)  to 
be  computed  at  a  simple  bifurcation  point,  with  one  solution  corresponding  to  each 
branch.  Using  the  appropriate  pair  of  (lio,  Ao)  in  (2.8)  allows  brand)  switching.  In  [7] 
a  more  detailed  study  of  the  singular  behavior  and  brand)  switching  ax  bifurcation  is 
given. 

In  order  to  solve  the  linear  system  in  (2.10)  by  direct  methods,  several  approaches 
are  possible.  One  way  is  to  perform  Gaussian  elimination  on  the  inflated  matrix  A, 
with  some  form  of  pivoting  to  ensure  stability.  But  this  approach  completely  ignores 
the  sparse  structure  which  is  usually  found  in  G.  arising  from  nonlinear  elliptic 
eigenvalue  problems.  In  order  to  take  advantage  of  the  structure  in  Cm  Keller  [12] 
suggested  the  following  block -elimination  procedure: 

Algorithm  BE  (block -elimination) 

Solve 

(2.11)  G„y  ■  G* 
and 

(2.12)  G^--G. 

Set 

(2.13)  SA  -  (Slz-S)/(S\  -  Sly) 

and 

(2.14)  6u-:-6A  y. 
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Note  that  only  systems  with  the  coefficient  matrix  G„  have  to  be  solved,  so 
structures  in  Cu  can  be  exploited.  Moreover,  only  one  factorization  of  G_  is  needed. 
It  has  been  shown  [27]  that  even  when  G„  is  becoming  singular.  Algorithm  BE 
produces  iterates  that  converge  quadratically  at  limit  points. 

Continuation  methods  of  various  forms  and  levels  of  sophistication  have  been 
widely  used  in  the  engineering  literature.  For  a  recent  survey  of  numerical  methods 
for  bifurcation  problems,  see  for  example  [IS].  The  approach  taken  here  is  due  to 
Keller  [12],  and  has  recently  been  applied  to  other  problems  in  fluid  mechanics  [5], 
[6],  [15],  [16],  [25],  [27].  A  related  approach  suggested  by  Abbott  [1]  corresponds 
(in  a  loose  way)  to  applying  Algorithm  BE  to  the  matrix  A  with  the  last  column 
permuted  into  the  first  n  columns  so  that  the  corresponding  coefficient  matrix  in 
equations  (2.11)  and  (2.12)  becomes  nonsingular  even  at  limit  points.  However,  as 
has  already  been  pointed  out,  any  structure  or  symmetry  in  G„  is  lost  in  the  process, 
and  hence  that  approach  seems  unsuitable  for  large  elliptic  systems  in  two  or  three 
dimensions. 

3.  Multi-grid  methods. 

3.1.  introduction.  The  class  of  multi-grid  (MG)  methods  that  we  use  here  is 
based  on  work  by  Bakhvalov  [3],  Brandt  [4],  Federenko  [8],  Hackbush  [10]  and 
Nicolaides  [19],  We  shall  only  briefly  describe  here  the  particular  MG  algorithms  that 
we  have  used  for  linear  elliptic  problems  that  arise  in  our  treatment  of  nonlinear 
elliptic  eigenvalue  problems. 

The  particular  way  in  which  we  use  the  MG  idea  is  to  use  a  hierarchy  of  grids, 
rather  than  a  single  one,  in  order  to  speed  up  the  convergence  rate  of  the  solution 
process.  The  MG  process  has  some  very  desirable  theoretical  properties:  for  certain 
elliptic  operators  on  an  n  x  n  grid,  it  computes  the  approximate  solution  to  truncation 
error  accuracy  in  0(n3)  arithmetic  operations  and  0(n3)  storage.  It  seems  natural  to 
consider  the  use  of  MG  methods  for  solving  nonlinear  eigenvalue  problems.  MG 
methods  have  been  applied  to  solution  of  linear  eigenvalue  problems  by  Hackbush 
[11]  and  McCormick  [17]. 

3.2.  The  Cyde  C  V  ^  algorithm.  The  particular  MG  algorithm  that  has  been 
used  in  this  study  is  based  or,  the  "Cycle  C"  algorithm  described  in  Brandt  [4],  This 
is  an  algorithm  for  iteratively  solving  the  discrete  equations  approxiir-ting  a  linear 
elliptic  problem  on  a  given  grid,  through  interaction  with  a  hierarchy  of  coarser  grids, 
taking  advantage  of  the  fact  that  the  different  discretisations  on  the  different  grids 
are  all  approximations  to  the  same  continuous  problem.  We  note  that  there  are  other 
MG  algorithms  [4]  proposed  for  implementing  continuation  procedures  outside  of  the 
context  of  the  pseudo  arc-length  framework.  Some  potential  problems  with  these 
related  algorithms  are  discussed  in  §  3 .4.  We  do  not  know  how  well  such  MG  algorithms 
perform  and  we  hope  to  carry  out  our  own  investigation  on  such  related  methods  in 
the  future.  In  this  paper,  MG  algorithms  are  used  to  solve  the  fine  grid  discrete 
equations  that  arise  in  the  pseudo  arc-length  continuation  procedure. 

Consider  a  hierarchy  of  grids  (G°,  C\  •  •  • ,  GM),  with  GM  being  the  finest  one, 
defined  on  a  domain  ft  with  corresponding  mesh  sizes  ( A0  >  A i  >  •  •  •>/»*,).  and  all 
approximating  the  same  linear  elliptic  problem: 

(3.1)  LU  *F  on  ft,  Um0  on  5ft. 

The  discrete  equation  on  a  grid  G‘  is  written  as: 

(3.2)  LkUk-Fk  on  G\  =  0  on  5ft. 
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We  are  primarily  interested  in  obtaining  the  approximating  solution  Uu  on  the  finest 
grid,  and  we  shall  start  with  an  initial  guess  on  G"  and  apply  a  standard  relaxation 
procedure  such  as  the  Gauss-Seidel  procedure.  It  is  well  known  that  the  error  is 
reduced  rapidly  in  the  first  few  iterations  but  then  the  reduction  rate  becomes  very 
slow.  By  a  frequency  analysis,  it  can  be  shown  that  the  fast  reduction  occurs  when 
the  residual  (or  the  error)  in  the  current  iterate  has  large  harmonics  on  the  scale  of 
the  grid,  the  so-called  high  frequencies.  Now  at  a  stage  in  the  iterative  process  where 
the  error  reduction  rate  slows  down,  let  the  current  iterate  be  u*.  Define  the  error 
v*  in  the  iterate  as  v*4  ■  UM Then  the  error  v"  satisfies  the  following  equation: 

(3.3)  L“t“-F“-L“u“-R“  on  G*.  e"«0  on  *GM. 

The  residual  R*  is  computable,  and  hence  the  original  problem  of  solving  for  UM 
can  be  reduced  to  an  equivalent  one  of  solving  (3.3)  for  t>M.  There  seems  to  be  no 
obvious  advantage  in  using  (3.3)  rather  than  continuing  with  the  original  relaxation 
procedure  with  uM.  However,  if  the  error  t?M  and  the  residual  R*  are  smooth  relative 
to  GM,  that  is,  if  their  high  frequency  components  have  been  smoothed  out  by  the 
relaxation  procedure,  then  we  can  approximate  the  solution  of  (3.3)  on  a  coarser  grid. 
My  GM~\  by  solving: 

on  G“-\ 

(3.4) 

v  *0  on  oG  , 

After  this  problem  is  solved  we  can  interpolate  the  solution  vu~l  onto  G*  to  get: 

(3.5)  new  u*  -  old  +  w*-i  /JJ- 1***'1, 

where  wm-i  is  an  interpolation  factor,  normally  taking  the  value  unity,  and  I\  sunds 
for  some  interpolation  operator  from  G‘  to  G'.  The  solution  process  for  equation 
(3.4)  on  Gm~'  usually  costs  considerably  less  than  the  cost  of  solving  equation  (3.3) 
on  GM.  If  v*  is  indeed  smooth  (relative  to  GM),  then  GM_1  should  proride  adequate 
resolution  for  v *  and  hence  iti- \vM~l  should  be  a  good  approximation  for  t  M.  This 
principle  of  transferring  to  a  coarser  grid  when  convergence  slows  down  can  be  applied 
recursively.  Thus  for  example,  we  can  start  with  a  zero  initial  guess  for  i>M_1  in  equation 
(3.4)  and  apply  the  Gauss-Seidel  relaxation  procedure  to  the  iterates  on  GM-1.  When 
the  convergence  slows  down,  we  can  again  transfer  to  the  next  coarser  grid  GM'J, 
and  so  on.  One  can  view  the  whole  process  as  each  grid  smoothing  just  those  frequencies 
in  the  error  that  are  high  relative  to  its  own  mesh  size,  each  doing  its  job  efficiently 
because  these  high  frequencies  are  precisely  those  that  are  efficiently  smoothed  out 
by  relaxation  procedures. 

The  control  of  when  to  transfer  between  grids  can  follow  a  fixed  strategy  or  an 
adaptive  one.  A  fixed  strategy  could  be  of  the  following  kind  (see  Nicolaides  [19]): 
perform  p  relaxation  sweeps  on  each  grid  Gk  before  transferring  to  a  coarser  grid 
G*~\  and  perform  q  relaxation  sweeps  before  interpolating  back  to  a  finer  grid  G**1. 
An  adaptive  strategy  could  be  as  follows  (tee  Brandt  [4]):  transfer  to  a  coarser  grid 
when  the  ratio  of  the  residual  norm  of  current  iterate  to  the  residual  norm  a  sweep 
earlier  is  greater  than  some  tolerance  tj,  and  transfer  to  a  finer  grid  when  the  ratio 
of  the  residual  norm  of  current  iterate  to  the  residual  norm  on  the  aext  finer  grid  is 
leu  than  another  tolerance  6.  For  simple  problems  like  Poisson’s  equation  on  a  square, 
the  overall  MG  efficiency  is  very  insensitive  to  which  particular  strategy  is  used  and 


180 


TONY  F.  C.  CHAN  AND  H.  B  KELLER 


what  values  are  used  for  (p,  q)  or  (vj,  8).  We  shall  refer  to  the  above  particular  fixed 
strategy  the  ( p ,  q)  strategy  and  the  adaptive  strategy  the  (tj,  6)  strategy. 

3.3.  Indefinite  problems.  In  the  Cycle  C  algorithm  just  described,  convergence 
on  the  lowest  (coarsest)  grid  G°  is  obtained  by  repeated  relaxation  sweeps.  For  positive 
definite  matrices,  convergence  on  G°  can  be  guaranteed.  For  indefinite  problems, 
however,  convergence  on  G°  cannot  be  obtained  by  repeated  relaxation  sweeps, 
because  the  components  of  the  error  that  correspond  to  eigenfunctions  with  negative 
eigenvalues  will  grow  as  a  result  of  relaxation  sweeps  (see  the  analysis  in  §  5).  Therefore, 
for  indefinite  problems,  a  direct  solution  (e.g.,  Gaussian  elimination)  must  be  employed 
on  the  coarsest  grid.  If  this  coarsest  grid  is  fine  enough,  it  will  also  provide  corrections 
to  those  growing  components  of  the  iterates  on  all  finer  grids.  However,  too  fine  a 
grid  for  G°  will  increase  the  cost  of  the  direct  solution  procedure.  Hence  a  little  care 
must  be  taken  regarding  the  size  of  the  coarsest  grid  for  indefinite  problems.  Fortu¬ 
nately,  for  “not  too  indefinite"  problems  G°  can  be  chosen  coarse  enough  so  that  the 
direct  solution  on  G°  will  not  aflect  the  overall  efficiency  of  the  MG  procedure 
seriously.  Since  indefinite  problems  occur  frequently  in  nonlinear  elliptic  eigenvalue 
problems  and,  in  particular,  in  our  model  problem,  we  shall  use  such  a  direct  solution 
on  G°  whenever  necessary. 

3.4.  Continuation  methods.  Brandt  [4]  suggested  using  continuation  methods  in 
conjunction  with  the  MG  procedure.  His  main  idea  is  to  use  coarse  grids  for  continu¬ 
ation,  with  little  work  and  crude  accuracy,  and  only  use  the  finer  grids  at  the  final 
continuation  step  to  achieve  higher  accuracy.  We  have  not  pursued  this  idea  here. 
We  believe  that  it  will  work  as  long  as  we  stay  away  from  singular  points.  Around  a 
limit  point,  however,  the  solution  branches  corresponding  to  different  grids  may  look 
like  the  situation  in  Fig.  3.1.  If  we  continue  on  the  coarse  grid  to  A  *  and  try  to  refine 


using  the  finer  grid,  while  keeping  A*  fixed,  we  cannot  hope  to  obtain  a  fine  grid 
solution  because  A*  is  larger  than  the  fine  grid  limit  point  Af  (i.e.,  no  fine  grid  solution 
exists  for  A  >A/).  In  the  opposite  case,  there  is  no  coarse  grid  solution  at  A*  so  we 
cannot  get  started  on  that  grid.  Hence,  in  general,  we  have  to  be  extremely  careful 
in  using  MG  methods  and  continuation  around  singular  points. 
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4.2.  Arc-length  continuation  with  direct  methods.  We  first  apply  the  arc-length 
continuation  method  of  §  2  to  (4. 1 )  using  direct  methods.  For  this  problem,  a  trivial 
solution  is  lu  *0,A» 0).  We  can  thus  start  at  this  trivial  solution  on  the  lower  branch 
and  march  along  the  solution  branch,  past  the  limit  point,  and  continue  on  to  the 
upper  solution  branch.  Since  the  only  singular  point  in  this  problem  is  a  limit  point, 
this  in  principle  presents  no  problem  to  the  arc-length  continuation  procedure, 
although  the  step  size  might  have  to  be  reduced  and  controlled  appropriately  near 
the  limit  point.  If  desired,  the  limit  point  can  be  accurately  determined  by  other  related 
techniques  [1],  [13]. 

The  derivatives  of  the  operator  G  in  equation  (4.1)  that  are  needed  for  the 
arc-length  continuation  technique  are: 

(4.2)  Gu  =  A  + A  e“, 

(4.3)  G*  ■  e“. 

Now  if  we  approximate  the  Laplacian  operator  by  the  standard  five-point  stencil  on 
a  uniform  grid,  the  operator  Cu  will  be  approximated  by  the  usual  block  tridiagonal 
matrix  and  the  operator  G*  by  a  column  vector. 

In  the  application  of  the  arc-length  continuation  technique,  we  will  have  to 
repeatedly  solve  linear  systems  of  equations  with  the  matrix  given  by  G„.  The  solution 
of  these  linear  systems  is  the  central  part  of  the  arc-length  continuation  method. 
Hence,  an  efficient  linear  system  solver  is  crucial  to  the  overall  performance  of  the 
continuation  technique.  In  this  section,  we  present  some  computational  results  for 
Bratu's  problem  using  a  direct  method  (Gaussian  elimination)  of  solution  of  the 
linearized  difference  equations.  For  large  problems,  this  would  be  prohibitively  expen¬ 
sive.  However,  the  results  here  are  intended  to  demonstrate  the  performance  of  the 
continuation  procedure  independent  of  the  linear  algebra  method  employed.  In  the 
next  section,  we  shall  investigate  the  use  of  multi-grid  methods  for  solving  the  linear 
equations.  It  should  be  pointed  out  that  G„  is  generally  not  separable,  and  therefore 
we  cannot  use  fast  Poisson  solvers  directly  even  on  rectangular  domains.  Moreover, 
this  matrix  is  indefinite  on  the  upper  branch,  and  hence  iterative  methods  like 
successive-over-relaxation  cannot  be  used  directly. 

We  present  some  of  our  computed  results  in  Table  4.1  and  Fig.  4.2.  Only  the 
behavior  of  the  solution  branch  near  the  limit  point  for  a  few  relatively  coarse 
discretizations  is  presented.  This  is  to  be  compared  with  the  values:  A*  * 6.8081 1698 
and  u(.5,  .5)=“  1.3916603  for  a  grid  with  h  *  j1?  with  the  nine-point  finite  difference 
operator  as  computed  by  Abbott  [1]  and  to  the  easily  obtainable  exact  solution 
(A*  ■  18/e  *6.62183,  u*  *  1)  for  the  case  h  =  As  expected,  the  step  size  ds  *=s-sr, 
had  to  be  suitably  controlled  near  the  limit  point,  but  otherwise  we  encountered  no 
difficulty  in  continuing  past  the  limit  point. 

4  J.  Arc-length  continuation  with  multi-grid  methods.  In  this  section  we  discuss 
the  use  of  MG  methods,  rather  than  direct  methods,  for  solving  the  linear  equations 
that  arise  in  the  continuation  procedure.  The  MG  method  that  we  use  was  described 
in  §  3  and  Gauss-Seidel  is  the  smoothing  relaxation  process.  Since  the  Jacobian  matrix 
Gu  becomes  indefinite  on  the  upper  branch,  we  use  a  direct  method  on  the  coarsest 
grid  in  the  neighborhood  of  the  limit  point  and  on  the  upper  branch. 

We  started  the  continuation  procedure  with  the  trivial  solution  (u  *  0,  A  *  0), 
with  h  ■  i  on  the  coarsest  grid,  and  a  total  of  four  levels  of  grids,  making  the  finest 
grid  with  h  -  As  expected,  the  MG  method  worked  fine  and  we  were  able  to 
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Table  4.1. 

Computed  results  for  Brotu's  problem  near  limit  point 


6.000000 

0.619061 

6.485170 

0.809435 

6.572858 

0.883052 

6.621830 

0.999899 

6.614022 

1.04937 

.8885E -4  —  limit  point 


6.500000 

1.00456 

0.9632 

6.689007 

1.14350 

0.9041 

6.802681 

1.34995 

0.2965 

6.805499 

1.39043 

-1.1732E-4  —limn  point 

6.805485 

1.39368 

-0.0125 

FlG.  4.2.  Computed  results  for  Bratu  s  problem  near  limit  point. 


continue  up  to  very  close  to  the  limit  point,  at  A  *  6.604  on  the  lower  branch.  However, 
we  noticed  that  the  convergence  of  the  MG  method  deteriorates  as  we  move  in  towards 
the  limit  point.  For  example,  the  number  of  equivalent  relaxation  sweeps  on  the  finest 
grid  required  to  reduce  the  residua)  norm  by  an  order  of  magnitude,  which  is  a 
convenient  way  of  measuring  the  efficiency  of  MG  methods,  went  from  about  5  at 
A  •  0  to  about  20  at  A  ■  6.803  and  to  divergence  at  A  *  6.805.  The  divergence  occurred 
in  the  MG  method  and  not  in  the  Newton  iteration.  It  is  not  due  to  the  possible 
indefiniteness  of  the  Jacobian  matrix  on  the  finest  grid.  This  can  occur  near  the  limit 
point  after  a  large  Euler-predictor  step.  We  performed  other  tests  starting  on  the 
upper  branch,  away  from  the  limit  point,  where  the  Jacobian  matrix  is  indefinite,  and 
here  the  MG  method  performed  as  efficiently  as  on  the  lower  branch.  From  our 
experience,  this  divergence  is  strictly  a  phenomenon  associated  with  the  limit  point. 


w 


3C 
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and  to  the  best  of  our  knowledge,  has  never  been  discussed  or  analyzed  in  the  literature. 
We  study  this  effect  in  §  5. 

The  exact  value  of  A  at  which  this  divergence  first  occurs  varies  slightly  with  the 
size  of  the  coarsest  grid  h0,  but  is  quite  independent  of  the  other  parameters  of  the 
Cycle  C  algorithm  (e.g.,  tj  and  S).  In  all  the  cases  we  have  run,  this  divergence  made 
it  impossible  to  continue  past  the  limit  point.  Therefore,  a  remedy  is  needed.  Before 
we  can  find  one,  we  must  understand  the  reason  for  the  divergence. 

5.  Analysis  of  multi-grid  methods  for  near-singular  systems.  For  the  present 
analysis,  we  assume  that  the  linear  operator  L  is  self-adjoint  and  has  the  complete 
set  of  orthonormal  eigenfunctions  •  •  •}  with  corresponding  real  eigenvalues 

=  "  }•  The  operator  G„  in  the  Bratu  problem  clearly  satisfies  the  above 

hypothesis.  Thus  the  solution  U  to  LU  =  F  can  be  written  as: 


(5.11 


U=Za,£„  /  =  1,2,' 


We  assume  that  the  discrete  approximations  Lk  to  the  continuous  L  are  symmetric. 
Thus  they  have  real  eigenvalues  f  S  m  :  S  •  •  •  S  m  v* }  and  a  complete  set  of  orthonor¬ 
mal  eigenvectors  {f  {,  £;,••■,  £s, }.  Here  Nk  is  the  dimension  of  the  matrix  representing 
Lk.  For  most  reasonable  approximations,  and  certainly  for  the  five-point  formula  used 
for  the  Bratu  problem  on  a  rectangle,  this  is  true. 

Assume  that  after  iterating  (relaxing)  on  the  grid  Gk,  convergence  has  slowed 
down  and  a  transfer  to  the  next  coarser  grid  is  desired.  Let  the  current  iterate  be  uk, 
and  the  corresponding  "correction"  be  t*  so  that  Uk  *  uk  + vk  where  Uk  satisfies 
LkUk  *  Fk.  The  correction  problem  is  given  (as  in  §  3)  by: 


if.  2) 


Lkvk  =  Rk~Fk-Lkuk  in  Gk, 


ck  =0  on  dGk. 


This  is  approximated  on  G  by 

(5 .3 1  Lk'\k~'  =  /£-,**  in  Gk,  r*''  =  0  on  dGk~l. 


Using  the  eigenvector  expansion  of  v  in  (5.2)  we  get 

v. 


(5.4) 

where 

(5.5) 


a , 


v 


*  <R\  f*> 


V  nk-k 

1  a.c.. 


i  -  1.  •  ■  ■ ,  A>*. 


Suppose  now  that  (5.3)  is  solved  exactly  (by  either  direct  solution  or  Cycle  C  or  any 
other  means)  on  Gk~\  The  solution  vk~'  is  then 


(5.6) 

where 

(5.7) 


.‘-i  =  \ 


I  ak-'(k-\ 

i- 1 


{Ikk-'Rk.(k-') 
a.  - - n - • 

Ml 


MULTI-GRID  CONTINUATION 


Tbe  key  idea  in  the  MG  method  is  that  if  vk  and  Rk  are  smooth  enough,  they  can  be 
well  approximated  on  Gk~\  Thus  it  is  important  for  efficiency  considerations  that1 

(5.8)  Ikk-ivk"mvk. 

% 

Using  (5.4)  and  (5.6),  this  is  equivalent  to: 

(5.9)  *1  V-'/J-.f*- ‘  »  I  aUl 


be  the  case  if 

*4-1 

2.  o,  /*- if, 

l*’  *  i 
«- 1 

(a) 

lfi/SN* 

(b) 

flf'1  »fl,\ 

1S/SN* 

(c) 

(5.10)  (a)  /I- iff  "*?.  lS/Stf*-,, 

(5.11)  (b)  of'1  mak,  1  fiSATfi. 

(5.12)  (c)  aim  0,  />**-.. 

Conditions  (5.10)  and  (5.11)  ensure  that  the  coarse  grid  correction  vk~l  improves  the 
lower  modes  of  the  iterate  uk.  Condition  (5.12)  is  essentially  the  smoothness  required 
of  vk  on  Ck  (i.e.,  negligible  higher  modes). 

Now  condition  (5.10)  is  satisfied  for  the  low  frequency  eigenfunctions  of  the 
continuous  operator  L  if  the  grids  Ck  and  Gk~l  are  both  fine  enough  to  resolve  these 
eigenfunctions.  This  holds  in  many  cases  since  the  lower  eigenfunctions  of  most 
second-order  elliptic  operators  over  smooth  domains  are  very  smooth.  For  the  Bratu 
problem,  the  eigenfunctions  are  very  close  to  products  of  sines  and  cosines  (the 
eigenfunctions  of  the  Laplacian  operator)  and  so  the  lower  modes  arc  easily  resolved 
by  very  coarse  pids.  Condition  (5.11),  on  the  other  hand,  turns  out  to  be  violated  if 
the  operator  L*  is  near  singular.  This  is  what  caused  the  divergence  of  the  Cycle  C 
algorithm  in  the  arc-length  continuation  procedure  as  we  approach  the  limit  point 
(see  §  4.3).  We  shall  analyze  this  case  next. 

From  (5.5)  and  (5.7),  condition  (5.11)  becomes 


(5.13) 


- m - * — i — , 

M. 


iS/SN*-,. 


We  claim  that  if  condition  (5.10)  is  satisfied,  and  if  the  transfer  from  Ck  to  Gk~l  is 
done  only  after  the  residual  Rk  has  been  smoothed,  then  the  numerators  in  (5.13) 
will  have  approximately  the  same  value.  To  show  this,  we  expand  Rk  as 


(5.14) 

where 

(5.15) 


r  .*?. 


Thus  the  numerator  on  the  right-hand  side  of  (5.13)  is  precisely  r,.  To  estimate  the 
numerator  on  the  left  hand  side  of  (5.13),  we  proceed  as  follows: 

.  ,  **  ...  *4-1  ...  *4  ... 


(5.16) 


r/J-tf-  I  rjkk-'€i 


l  r/i'V.*. 


1  We  shall  me  the  ■  symbol  to  mean  rather  loosely  "approximately  equal  to"  The  meaning  should 
be  clear  by  context  Also,  we  shall  assume  that  the  interpolation  (actor  in  equation  »3.5t  is  equal  to 
one  unless  sated  otherwise. 
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Now  if  condition  (5.10)  holds,  its  converse 


(5.17)  1 S  iSN't-i 

also  holds.  Also,  if  has  been  smoothed  on  G\  then  r,  [for  .V*.,  <  /S.\*]  must  be 
small  compared  with  r,  [for  15  /  S\k-i].  Alternatively  (5.12)  assumes  a*  •  r,ink,  *0 
for  i>Sk- j.  Therefore,  we  can  approximate  in  (5.16)  by  dropping  the  second  sum 
on  the  right-hand  side  to  get 


(5.18) 


Hence 

(5.19) 


i-i 


(5.19)  </i'lJ?*,ff'*>*rB  IS  /SAT*.,. 

Therefore,  from  (5.15 1  and  (5.19).  we  have,  as  claimed  earlier. 


(5.20) 


[»  foris/sN*.,. 

The  relations  in  (5.20)  imply  that  condition  (5.13)  will  be  true  if 


(5.21) 


FT*1,  lSiSN*-,. 


Actually,  these  conditions  need  to  be  strengthened  in  order  to  guarantee  that  the  visit 
to  G‘~‘  actually  improves  the  accuracy  of  u*.  This  can  be  seen  as  follows.  The  error 
in  the  iterate  uk  before  the  transfer  to  G*-1  is  given  by 


(5.22) 


old  error*  vh  *  I  a *£* . 


From  (3.5),  the  new  error  in  uk  after  coming  back  from  a  visit  to  Gk~l  is  given  by 
(5.23)  new  error*  vk -*•*- 

In  view  of  (5.41  and  (5.6),  the  above  gives 


(5.24) 


N._, 

new  error  a  £  (a  *  -  wk  .  \a  *  “ 1 ){ ,k  +■  higher  modes 

c-l 

a  X  ( 1  — — )  a  f  +  higher  modes. 

i-l  \  a*  / 


From  (5.5),  (5.7)  and  (5.20),  we  have 


and  therefore  we  can  write  the  new  error  in  (5.24)  as 

(5.25)  new  error  a  £  ’  ( 1  -  —-y'.T1  )  a  ff  f  +  higher  modes. 

For  obvious  efficiency  and  convergence  considerations,  the  new  error  should  preferably 
be  less  than  the  old  error,  at  least  for  the  lower  modes.  In  other  words,  condition 
(5.21)  should  be  strengthened  to 


(5.26) 


1— <1. 
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i.e., 

(5.27)  *  0<— for  lSiSN*.,. 

Now  if  the  ratios  of  eigenvalues  in  (5.21)  are  not  dose  to  unity,  the  interpolation 
factors,  *•>*_,,  should  be  chosen  so  that  condition  (5.27)  is  satisfied.  Otherwise  the 
new  error  can  be  larger  than  the  old  error  in  some  modes. 

It  should  be  pointed  out  that,  in  general,  condition  (5.27)  is  not  necessary  tor 
the  convergence  of  the  “Cycle  C“  algorithm.  This  is  the  case,  for  instance,  if  L  and 
the  L*’s  are  all  positive  definite.  Then  Gauss-Seidel  sweeps  on  any  grid  Gk  will  reduce 
the  amplitude  of  every  mode  present  in  the  error.  In  such  cases,  convergence  on  any 
grid  can  be  achieved  by  merely  doing  enough  relaxation  sweeps.  Then  it  is  not  necessary 
for  the  next  coarser  grid  to  provide  any  improvement  on  the  current  iterate,  although 
it  would  obviously  improve  the  efficiency  of  the  overall  algorithm  if  it  does  so.  In  fact, 
the  MG  method  derives  its  effidency  from  the  very  fact  that  the  coarser  grids  do 
provide  improvements  in  the  current  iterate  uk  in  the  lower  modes.  These  are  precisely 
those  modes  that  have  poor  convergence  rates  for  the  relaxation  sweeps  on  <7*.  Thus, 
even  in  the  positive  definite  case,  it  is  important  (from  an  efficiency  viewpoint)  that 
conditions  (5.27)  hold,  at  least  for  small  /’ s. 

If  the  operator  L  and  the  Lk't  are  indefinite  the  situation  is  different  because 
some  modes  will  grow  if  we  simply  perform  relaxation  sweeps  on  a  fixed  grid.  Such 
modes  have  to  be  corrected  by  going  to  coarser  grids  and  using  a  direct  method  on 
the  coarsest  grid.  Further,  the  interpolation  factors,  *»■*_,,  should  be  chosen  such  that 
condition  (5.27)  is  satisfied  for  these  modes.  Condition  (5.27)  has  been  suggested  by 
Brandt  [4]  for  indefinite  problems.  However,  as  we  show  later,  most  nonlinear 
eigenvalue  problems  with  limit  points  and  bifurcation  points  abound  with  indefinite 
operators,  but  they  do  not  cause  difficulties  in  the  sense  of  violating  condition  (5.27). 
Essentially  only  one  mode  causes  problems  on  each  Ck  and  it  is  the  mode  that 
corresponds  to  the  eigenvalue  that  is  nearest  zero  as  the  singular  point  is  approached. 
Merely  including  the  interpolation  factors  so  that  condition  (5.27)  is  satisfied  turns 
out  to  be  very  inefficient.  Further,  it  is  not  clear  that  such  factors,  can  be  found 
at  all  in  this  case. 

Another  source  of  difficulty  is  that  the  process  of  interpolating  vk~i  into  Ck 
introduces  high  frequency  errors.  That  is,  the  exact  relation  corresponding  to  (5.10)  is: 

(5.28)  +  «-l, 2, •••,**.„  forlSiSAVi. 

/-i 

and  the  coefficients  bk  may  be  large  for  This  would  result  in  a  violation  of 

(5.12).  Fortunately,  these  high  frequency  errors  are  very  efficiently  smoothed  out  by 
the  subsequent  relaxation  sweeps  on  Gk,  and  thus  these  errors  are  automatically 
corrected. 

For  elliptic  operators  which  are  “far"  from  being  singular  and  with  a  reasonable 
grid  system  {<?*}  condition  (5.27)  can  be  assured.  For  example,  if  Z.  is  the  negative 
Lapladan,  —A,  on  a  unit  square  with  Dirichlet  boundary  conditions,  then  it  is  known 
(e.g.,  [9])  that  the  eigenvalues  of  L  are  given  by 

(5.29)  Mm.»  “  (mn)7  (nirf. 

The  corresponding  eigenfunctions  are: 

(5.30)  €m.n  "Sin  (mm)  sin  (str y). 


These  eigenfunctions  evaluated  at  the  discrete  interior  grid  points  of  a  uniform  mesh 
on  the  unit  square  give  the  eigenfunctions  of  the  discrete  S -point  approximations, 
Lk  *  -A*,  with  h  being  the  uniform  mesh  size.  The  eigenvalues  of  Lk  are,  with 
Sx  «5y  ■  ht, 

/r  ,  k  ^[sin5  (m»rfik/2)+sinJ  (nirhk/2)] 

(5.31)  ■ - n - • 

fit 

Some  of  these  eigenvalues  are  tabulated  in  Table  5.1  for  various  mesh  sizes,  hk.  The 
ratios  are  given  in  Table  5.2.  We  see  from  Table  5.2  that  condition  (5.27) 


Table  S.l. 
for  -A* 

km 

0 

1 

2 

3 

ao 

(m,n) 

hom\ 

maM 

h,-i 

A«-0 

1.1 

16.0 

18.743 

19.487 

19.676 

19.739 

2.1 

NA 

41.37258 

47.238 

48.812 

49.348 

1.2 

NA 

41.37258 

47.238 

48.812 

49.348 

2.2 

NA 

64.0 

74.981 

77.947 

78.957 

3.1 

NA 

NA 

88.760 

96.126 

98.696 

1.3 

NA 

NA 

88.760 

96.126 

98.696 

3.2 

NA 

NA 

116.507 

125.261 

128.305 

2.3 

NA 

NA 

116.507 

125.261 

128.305 

3,3 

NA 

NA 

158.033 

172.575 

177.653 

Table  52. 

Ratios  for  -A* 

(m.  n) 

1 

• 

1,1 

1.17 

1.04 

1.01 

2.1 

NA 

1.14 

1.03 

1,2 

NA 

1.14 

1.03 

2.2 

NA 

1.17 

1.04 

3.1 

NA 

NA 

1.08 

1,3 

NA 

NA 

1.08 

3.2 

NA 

NA 

1.08 

2.3 

NA 

NA 

1.08 

3.3 

NA 

NA 

1.09 

is  satisfied,  with  w*_i * 1,  for  all  lower  modes  shown.  These  ratios  are  very  close  to 
unity,  even  for  the  case  where  the  coarsest  grid  has  only  one  interior  point.  We  have 
seen  from  condition  (5.11)  that  this  closeness  to  unity  is  very  desirable  and  this  fact 
partly  explains  the  well-documented  success  of  MG  methods  for  the  Laplacian 
operator. 

Near  the  limit  point  of  the  Bratu  problem,  the  operator  L  ■  G»  *  A + Ac  “  behaves 
very  much  like  a  shifted  Laplacian  operator.  Gearly,  if  the  factor  e“  were  replaced 
fey  a  constant,  a  say,  then  Gu  is  replaced  by  the  Laplacian  operator  with  a  shift  aA. 
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Then  the  eigenvalue  ratio  n 
(5.32) 


i.i/mi.7*.  valid  for  aA  •  0.  is  replaced  by: 
Mi.l  ~aA 


Since  0  <u  <1.4,  the  factor  f*  does  not  vary  much  and  we  assume  this  approximation 
to  be  valid  for  some  a  >0.  The  situation  is  depicted  graphically  in  Fig.  3.1  for  the 
grid  system  that  was  used  for  Table  5.1.  As  the  shift  aA  approaches  the  group  of 
eigenvalues  corresponding  to  the  (1, 1 )  mode  from  below,  die  ratios  in  (5.31)  increase. 
As  aA  continues  to  increase,  the  ratio  of  eigenvalues  will  become  greater  than  2,  then 
increase  towards  +oc,  jump  to  -so  dis continuously,  and  start  increasing  from  -oo  to 
1.  The  situation  is  depicted  in  Fig.  5.2. 


aiiMdiut 


;  \ 
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*1-1-1.. I- 
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We  thus  see,  under  the  above  assumptions,  that  condition  (5.27)  is  first  violated 
by  the  lowest  mode  (i.e.,  the  (1, 1)  mode)  on  the  two  coarsest  grids  G°  and  G\  In 
fact  the  lowest  eigenvalues  for  the  Bratu  problem  computed  at  the  first  point  on  the 
solution  branch  where  Cycle  C  diverged,  yields  the  ratio  almost  exactly  2.  On  the 
other  hand,  even  at  this  point,  condition  (5.27)  is  satisfied  by  the  (1, 1)  modes  on  the 
finer  grids.  In  other  words,  the  divergence  of  Cycle  C  is  seen  to  be  caused  by  one 
near-singular  grid  out  of  the  whole  hierarchy  of  grids  present.  The  mode  that  becomes 
singular  at  the  limit  point  of  the  Bratu  problem  is  the  (1, 1)  mode,  and  this  occurs 
first  on  the  G°  grid.  As  the  limit  point  is  approached,  Lk  on  some  of  these  grids  may 
even  become  indefinite,  while  others  (the  finer  grids)  may  still  be  positive  definite. 
Essentially,  the  near-singular  grid  causes  the  (1, 1)  mode  component  of  the  correction 
vk~\  when  viewed  as  an  approximation  to  vk,  to  have  the  right  direction,  but  the 
wrong  magnitude.  This  phenomenon  is  not  limited  to  the  Bratu  problem.  The  only 
thing  special  about  this  problem  is  that  it  is  the  eigenvalue  of  the  (1, 1)  mode  that 
becomes  zero  at  the  limit  point.  For  other  problems,  the  eigenvalue  of  the  operator 
L  that  becomes  zero  as  the  singular  point  is  approached  might  correspond  to  other 
modes.  Although  the  singular  point  in  the  Bratu  problem  is  a  limit  point,  we  can 
expect  the  same  behavior  at  a  bifurcation  point. 

Having  now  understood  the  cause  of  the  divergence  of  the  MG  method,  in  the 
next  section  we  shall  discuss  some  modifications  to  the  basic  Cycle  C  algorithm  that 
are  designed  to  overcome  such  difficulties. 

6.  Remedies  aad  new  algorithms.  In  this  section  we  discuss  approaches  that  have 
been  devised  to  overcome  the  difficulties  with  the  MG  method  near  singular  points. 
The  first  goal  is  to  modify  the  basic  Cycle  C  algorithm  so  that  it  will  converge  for 
values  of  A  close  enough  to  the  limit  point  so  that  the  arc-length  continuation  procedure 
can  take  us  past  the  limit  point  onto  the  upper  solution  branch.  A  more  ambitious 
goal  is  to  modify  Cycle  C  further  so  that  it  will  converge  arbitrarily  close  to  the  singular 
point.  Such  an  algorithm,  when  used  in  conjunction  with  the  arc-length  continuation 
technique  for  tracing  solution  branches,  will  make  the  overall  algorithm  much  more 
robust.  Moreover,  such  an  algorithm  may  prove  to  be  useful  for  locating  singular 
points  accurately,  either  using  an  arc-length  continuation  based  procedure  [13],  or 
some  other  procedure  that  uses  the  operator  G„  near  the  singular  point  [22].  We  shall 
see  that  the  first  goal  is  relatively  easy  to  achieve,  whereas  the  second  goal  is  much 
more  difficult.  However,  we  have  devised  a  Cycle  C  based  algorithm  that  has  performed 
very  well  when  applied  very  close  to  the  limit  point.  The  approaches  that  we  have 
tried  and  that  lead  to  the  final  algorithm  will  be  discussed  in  this  section.  We  shall 
describe  them  in  the  sequence  that  they  were  tried. 

Before  we  proceed,  however,  we  have  to  explain  a  few  general  strategies  that 
were  used.  First  of  all,  Gauss-Seidel  and  many  other  relaxation  schemes  are  not  very 
effective  in  smoothing  the  lower  modes,  especially  modes  with  near-zero  eigenvalues. 
Hence,  these  modes  must  be  eliminated  by  means  other  than  relaxation,  even  on  the 
coarsest  grid.  Therefore,  unless  stated  otherwise,  we  shall  use  a  direct  solution  on  the 
coarsest  grid  even  though  the  operators  Lk 's  may  be  positive  definite.  This  does  not 
affect  the  overall  efficiency  very  much  because  the  coarsest  grid  has  so  few  points  that 
direct  solution  is  very  fast  and  efficient. 

Another  strategy  concerns  the  treatment  of  the  mode  that  causes  the  divergence, 
that  is,  the  mode  with  a  near -zero  eigenvalue,  say  £,.  In  all  the  algorithms  that  are 
discussed,  this  mode  is  treated  separately  from  the  other  modes.  To  do  this,  it  is 
essential  to  have  approximations  to  this  mode  and  to  its  corresponding  eigenvalues, 
say  it  and  du  respectively.  Here  we  have  to  strike  a  balance  between  accuracy  and 
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efficiency.  If  we  compute  the  exactly,  then  we  can  completely  eliminate  the  ft 
error  components  on  each  grid.  Thus,  the  problem  on  G*  can  be  reduced  to  one  in 
which  a{  is  zero  (see  (5.25)).  When  this  is  done,  we  do  not  need  to  satisfy  condition 
(5.27)  for  this  mode.  On  the  other  hand,  the  work  involved  in  computing  accurate 
approximations  to  #*t  and  ft  for  each  k  would  be  at  least  as  much  as  solving  the 
original  linear  system.  Our  compromise  has  been  to  compute  an  approximation  £?  to 
£ i  on  the  coarsest  grid,  G°,  by  a  few  steps  of  inverse  iteration  with  zero  shift  (since 
die  eigenvalue  we  want  is  near  zero).  This  it  very  inexpensive  since  Gc  is  quite  coarse 
and  the  LU  factors  of  L°  are  already  available.  Then  we  interpolate  ff  onto  the  finer 
grids.  To  eliminate  the  high  frequency  errors  introduced  in  these  interpolations,  we 
do  two  things:  (1)  use  higher  order  interpolation,  e.g.,  cubic  instead  of  linear;  (2) 
smooth  the  interpolated  eigenfunctions  by  performing  a  few  relaxation  sweeps  on 
L*f  t  *0.  Estimates  of  the  eigenvalues,  mi.  are  then  computed  using  the  Rayleigh 
quotients:  (ft,  L*f f ).  We  view  this  as  a  preprocessing  phase  of  the  algorithm  and  the 
extra  work  is  usually  minimal  compared  to  the  overall  work.  Furthermore,  since  the 
eigenfunctions  (not  the  eigenvalues)  do  not  change  very  much  in  the  neighborhood 
of  the  singular  points,  we  can  use  the  same  approximation  for  different  linearized 
operators  Lk.  The  storage  required  to  store  these  eigenfunctions  is  less  than  twice  the 
size  of  the  finest  grid. 

We  use  the  (n.d)  adaptive  version  of  the  Cycle  C  algorithm,  unless  otherwise 
stated.  The  first  modified  algorithm  is  the  following. 

6.1.  Under*  and  over-interpolation.  The  idea  is  to  choose  w*-i  in  (3.5)  for 
interpolation  onto  Gk,  such  that  condition  (5.27)  is  satisfied  for  f  t.  Clearly  the  value 


Mi 


is  in  some  sense  optimal  since  it  eliminates  the  fj  term  in  (5.25).  For  the  case  discussed 
in  (  4.3,  this  modification  allows  the  computation  to  continue  past  the  point  A  -  6.804, 
where  divergence  of  Cycle  C  first  occurred.  In  fact  (with  a  little  luck)  we  succeeded 
in  continuing  around  the  limit  point  onto  the  upper  branch.  Here  the  eigenfunction 
fj  no  longer  presented  difficulties  for  the  MG  algorithm.  For  some  of  these  cases  #t? 
is  actually  negative  and  therefore  (6.1)  yields  a  negative  value  for  W|.  In  this  case  the 
transfer  from  G°  to  G'  violates  condition  (5.27)  for  all  modes  other  than  £,.  The 
errors  in  these  modes  must  be  reduced  by  extra  relaxation  sweeps  on  Gl.  In  other 
words  G°  only  provides  a  proper  correction  on  G1  for  the  fi  mode,  all  higher  modes 
are  treated  incorrectly  during  the  transfer.  The  efficiency  of  the  algorithm  thus  suffers. 
This  effect  is  especially  pronounced  if  some  factors  w*  are  either  very  large  or  negative 
or  (worse)  both.  The  algorithm  is  very  sensitive  to  the  parameters  (rj.  6)  and  thus  is 
not  robust.  It  can  even  diverge  if  the  higher  modes  are  not  reduced  fast  enough  on 
G*  after  the  transfer  from  G  . 

Even  worse,  the  above  algorithm  will  not  work  for  indefinite  problems  in  which 
some  intermediate  eigenvalue  is  near  zero.  For  example,  if  the  spectra  of  the  Lk  are 
similar  to  those  in  Fig.  6.1,  the  interpolation  factors  wk  are  controlled  by  the  ft 
belonging  to  eigenvalues  mi  near  zero.  On  the  other  hand,  the  eigenfunctions  f‘i 
require  that  condition  (5.27)  be  satisfied  because  these  modes  cannot  be  liquidated 
by  relaxation.  Conflicts  can  occur  when  ft  requires  w*  to  be  negative  while  f‘-i 
requires  w*  to  be  positive.  Indefinite  problems  of  this  type  occur  frequently  in  nonlinear 
eigenvalue  problems.  Mere  under-  or  over-interpolation  must  run  into  difficulties  for 
such  problems,  near  the  singular  points. 
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I 

Origin 

FiC.  6.1.  Inurmtdiatt  eigenvalue  near  zero. 


The  above  considerations  make  it  clear  that  the  eigenfunction  with  the  near-zero 
eigenvalue  must  be  isolated  and  treated  differently  from  the  other  eigenfunctions.  We 
use  the  approximate  eigenfunctions  that  are  computed  in  the  preprocessing  phase  for 
this  purpose  in  the  following  procedure. 


6.2.  Under-  and  over-interpolating  the  singular  eigenfunction  only.  We  use  an 

interpolation  different  from  that  in  (3.5).  Specifically  if 


(6.2 1 


_ . 

*  1  o,  c. 

i-J 


on  Gk~\  we  interpolate  it  onto  Gk  by 

(6.3.  Yar’fr1. 


1-2 


Further  m*.i  is  chosen  to  satisfy  (6.1).  Since  we  only  have  an  approximation  to  ft, 
we  use,  instead  of  (6.3): 


<6.4. 


pk-i 


In  practice,  this  performed  much  better  than  indiscriminate  under-  and  over-interpola¬ 
tion  described  in  §6.1.  It  was  the  more  efficient  when  both  procedures  worked.  In 
many  cases  when  (6.1)  yields  large  and/or  negative  values  for  wk,  only  the  current 
scheme  converges.  In  principle,  it  will  also  work  for  indefinite  problems  like  that 
depicted  in  Fig.  6.1.  The  efficiency  in  most  cases  was  very  respectable:  in  the  range 
of  6-10  units  per  order  of  magnitude  reduction  in  the  residual.  It  is  also  quite  insensitive 
to  the  parameters  (» j,  6).  Thus,  it  can  be  used  very  efficiently  and  reliably  with  the 
arc-length  continuation  procedure  for  tracing  out  solution  branches. 

Unfortunately,  this  improved  algorithm  fails  when  the  magnitude  of  wk  becomes 
too  large.  This  occurs  when  Lk  is  very  nearly  singular,  that  is,  with  mi  very  close  to 
zero.  Since  we  only  have  an  approximation  f>  to  fj,  large  factors  wt  in  (6.4)  introduce 
very  large  errors  in  the  other  modes.  Moreover,  the  estimates  fik  using  Rayleigh 
quotients  tend  to  be  too  large  (relatively)  when  is  very  small.  Then  (6.1)  gives  a 
value  of  h*  that  is  too  small.  Both  of  the  above  result  in  lower  efficiency  and  reliability. 
In  extreme  cases,  this  makes  the  algorithm  impractical.  To  overcome  this  difficulty, 
we  devise  an  algorithm  that  will  work  even  if  one  of  the  operators  Lk  is  very  nearly 
singular.  For  this  we  employ  the  idea  of  skipping  a  grid. 


6.3.  Skipping  the  singular  grid.  The  previous  algorithm  fails  if  the  operator  is 
very  nearly  singular  on  one  of  the  grids,  say  Gk.  The  idea  here  is  to  simply  delete 
this  grid  from  the  hierarchy  of  grids  used  by  the  MG  algorithm.  If  the  remaining  grids 
are  not  as  singular  as  the  deleted  grid  it  would  seem  that  the  algorithm  described  in 
§  6.2  should  work.  However,  calculations  show  that  skipping  a  grid  can  cause  other 
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problems.  When  G‘  is  skipped,  the  mesh  changes  more  drastically  from  Gi_I  to 
G‘*\  and  hence  the  interpolation  in  (6.4)  (now  /Jl!  instead  of  /*-j)  introduces 
larger  errors  into  the  higher  modes  on  G**\  These  high-frequency  errors  can  cause 
divergence  of  the  MG  process  unless  controlled  properly  by  the  parameters  (q,  6).  A 
large  value  of  q,  say  between  0.6  and  0.9,  makes  the  algorithm  more  robust  but 
involves  more  work  than  for  a  smaller  value  of  %  say  0.5.  We  encountered  a  case 
where,  with  all  eise  the  same,  the  new  skipping  algorithm  converges  for  n  -  0.9  but 
diverges  for  rj  *  0.6.  Granted  with  i?  ■  0.9  the  algorithm  may  be  very  reliable,  such 
sensitivity  to  one  parameter  is  very  undesirable.  Therefore,  we  considered  the  following 
modification. 


6.4.  Skipping  the  singular  grid  for  the  singular  eigenfunction  only.  The  idea  is 
to  skip  the  singular  grid  G*  for  (t  only,  and  to  keep  it  for  smoothing  the  other  modes. 
In  the  actual  implementation,  we  modify  the  algorithm  described  in  §  6.2  to  use 


(6.5) 


w*_, 


for  (■.  and  w*_, « 1  for  all  other  modes  to  transfer  from  G*~'  to  Gk  and,  after  a  feu- 
smoothing  sweeps  on  G*,  transfer  to  G**1  with  w*  ■  1  for  all  modes.  Note  that  we 
do  not  try  to  solve  the  Gk  equations  for  t>\  Trying  to  do  that  would  result  in  large 
magnification  of  the  ff  component  in  vk,  since  is  near  aero.  This  would  in  turn 
cause  problems  during  the  transfer  to  G**1. 

In  addition,  we  have  experimented  with  using  a  mixture  of  the  adaptive  (n.  6) 
strategy  with  the  nonadaptive  (p,  q)  strategy  (cf.  13.2).  We  have  found  an  (t>.  <f) 
strategy  that  is  as  good  as  any  other  we  have  tried.  In  this  strategy,  we  use  tj  to  control 
when  we  terminate  relaxation  on  a  certain  grid  and  go  on  to  a  coarser  grid,  and  use 
q  to  control  how  many  sweeps  to  do  on  a  grid  after  transfer  from  a  coarser  grid  before 
interpolating  onto  a  finer  grid.  A  typical  set  of  parameters  that  worked  well  is  (tj »  0.6, 
q  «  2).  The  resulting  algorithm  is  fairly  insensitive  to  actual  values  of  tj  and  q  and  is 
quite  robust.  It  is  also  quite  efficient.  It  consistently  achieved  an  efficiency  of  less  than 
about  12  units  per  order  of  magnitude  reduction  in  the  residual  for  most  problems 
that  we  have  encountered.  Some  of  these  problems  have  very  singular  grids  which 
presented  difficulties  for  all  of  the  previous  algorithms. 


7.  Summary.  In  this  paper,  we  study  arc-length  continuation  techniques  and 
multi-grid  techniques  for  solving  nonlinear  elliptic  eigenvalue  problems.  We  have 
applied  these  techniques  to  solve  a  model  nonlinear  elliptic  eigenvalue  problem  (the 
Bratu  problem).  We  have  found  that  as  long  as  we  stay  away  from  singular  points, 
the  two  techniques  combined  to  give  a  very  powerful  and  efficient  procedure  for 
tracing  solution  branches.  Near  singular  points,  however,  the  standard  multi-grid 
method  has  difficulty  converging  on  the  linearized  elliptic  systems  that  arise  in  the 
continuation  procedure.  One  consequence  is  that  we  cannot  continue  past  the  limit 
point  in  the  model  problem.  This  divergence  is  successfully  analyzed  and  several 
modified  multi-gnd  algorithms  have  been  designed  based  on  this  analysis.  The  best 
of  these  modified  algorithms  performs  efficiently  and  reliably  arbitrarily  close  to  the 
singular  points.  This  enables  the  continuation  procedure  to  continue  past  the  limit 
point  with  no  difficulty.  It  seems  reasonable  that  this  modified  multi-grid  algorithm 
can  be  useful  in  more  general  situations  where  nearly  singular  elliptic  systems  arise, 
such  as  in  inverse  iteration  [11],  [17], 
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